# Compute growth ---------------------------------------------------------------
# Load 2008 census data and join additional species data
census_2008 <- bw_census_2008 %>%
left_join(bw_species, by = "sp")
# Load 2014 census data and filter out resprouts
census_2014 <- bw_census_2014 %>%
filter(!str_detect(codes, "R"))
bw_growth <-
# Merge both censuses and compute growth:
compute_growth(census_2008, census_2014, id = "treeID") %>%
# Convert data frame to sf object
st_as_sf(coords = c("gx", "gy")) %>%
# Clean variables
select(ID = treeID, species = trait_group, genus, dbh = dbh2, growth, geometry)
# From forestecology package vignette
# devtools::install_github("rvalavi/blockCV")
library(blockCV)
library(sfheaders)
scbi_2013 <- read.csv("data/scbi.stem2.csv") %>%
select(treeID, stemID, sp, quadrat, gx, gy, dbh, date = ExactDate, codes, status) %>%
mutate(date = lubridate::mdy(date)) %>%
filter(gx < 300, gy > 300, gy < 600)
scbi_2018 <- read.csv("data/scbi.stem3.csv") %>%
select(treeID, stemID, sp, quadrat, gx, gy, dbh, date = ExactDate, codes, status) %>%
mutate(date = lubridate::mdy(date),
dbh = as.numeric(dbh)) %>%
filter(gx < 300, gy > 300, gy < 600)
census_df1 <- scbi_2013
census_df2 <- scbi_2018
id <- "stemID"
scbi_growth_df <-
# Merge both censuses and compute growth:
compute_growth(census_df1, census_df2, id) %>%
# Convert growth from cm to mm to make result comperable
mutate(growth = growth/10,
sp = as.factor(sp))
# Add spatial information
cv_fold_size <- 100
max_dist <- 7.5
scbi_study_region <-
#tibble(x = c(0,400,400,0,0), y = c(0,0,640,640,0)) %>%
tibble(x = c(0,300,300,0,0), y = c(300,300,600,600,3)) %>%
sfc_polygon()
# Add buffer variable to data frame
scbi_growth_df <- scbi_growth_df %>%
add_buffer_variable(direction = "in", size = max_dist, region = scbi_study_region)
scbi_cv_grid <- spatialBlock(
speciesData = scbi_growth_df, theRange = 100, yOffset = 0.9999, k = 9, verbose = FALSE)
# Add foldID to data
scbi_growth_df <- scbi_growth_df %>%
mutate(
foldID = scbi_cv_grid$foldID
)
# Visualize grid
scbi_cv_grid$plots +
geom_sf(data = scbi_growth_df, aes(col=factor(foldID)), size = 0.1)
scbi_cv_grid_sf <- scbi_cv_grid$blocks %>%
st_as_sf()
focal_vs_comp_scbi <- scbi_growth_df %>%
create_focal_vs_comp(max_dist, cv_grid_sf = scbi_cv_grid_sf, id = "stemID")
metadata <- read_csv("data/scbi.spptable.csv")
focal_vs_comp_scbi_wide <- focal_vs_comp_scbi %>%
group_by(focal_ID, focal_sp, dbh, foldID, growth, comp_sp) %>%
summarize(comp_basal_area = sum(comp_basal_area)) %>%
pivot_wider(names_from = comp_sp, values_from = comp_basal_area, values_fill = 0)
head(focal_vs_comp_scbi_wide)
## # A tibble: 6 x 54
## # Groups: focal_ID, focal_sp, dbh, foldID, growth [6]
## focal_ID focal_sp dbh foldID growth caca ceoc fram frni havi litu
## <int> <fct> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 4 nysy 136. 4 0.103 1.31 0.849 0.204 1.44 0.642 28.3
## 2 5 havi 88 4 0.150 4.58 0.849 0 1.44 2.31 0
## 3 79 tiam 477. 3 -0.161 4.52 0 3.24 0 0 0
## 4 80 caca 51.5 4 0.253 2.45 0 0 0 0 13.8
## 5 96 libe 23 3 0.262 1.53 0 4.51 0 0 0
## 6 101 litu 654. 5 0.552 0 0 0 0 0 31.4
## # … with 43 more variables: nysy <dbl>, saal <dbl>, tiam <dbl>, ulam <dbl>,
## # cagl <dbl>, unk <dbl>, quru <dbl>, cato <dbl>, fagr <dbl>, libe <dbl>,
## # astr <dbl>, frpe <dbl>, quve <dbl>, ulru <dbl>, ceca <dbl>, caco <dbl>,
## # cofl <dbl>, qual <dbl>, juni <dbl>, acru <dbl>, caovl <dbl>, amar <dbl>,
## # prav <dbl>, vipr <dbl>, prse <dbl>, pist <dbl>, acne <dbl>, qupr <dbl>,
## # elum <dbl>, juci <dbl>, saca <dbl>, cade <dbl>, crsp <dbl>, ploc <dbl>,
## # chvi <dbl>, ilve <dbl>, divi <dbl>, crpr <dbl>, qumi <dbl>, qufa <dbl>,
## # quco <dbl>, rops <dbl>, pivi <dbl>
# write_csv(focal_vs_comp_scbi_wide, "data/focal_vs_comp_scbi_wide.csv")
Plot all trees in SCBI plot:
ggplot(data = scbi_growth_df) +
geom_sf(aes(col = sp), size = 0.5) +
labs(title = "All trees in SCBI site")
Plot all the trees (??)
ggplot(data = bw_growth %>% filter(species == "maple")) +
geom_sf(aes(col = species), size = 0.5)+
labs(title = "All maple trees in Michigan Big Woods site")
ggplot(data = bw_growth %>% filter(species == "evergreen")) +
geom_sf(aes(col = species), size = 0.5)+
labs(title = "All maple trees in Michigan Big Woods site")
ggplot(data = bw_growth %>% filter(species == "oak")) +
geom_sf(aes(col = species), size = 0.5)+
labs(title = "All maple trees in Michigan Big Woods site")
ggplot(data = bw_growth %>% filter(species == "short_tree")) +
geom_sf(aes(col = species), size = 0.5)+
labs(title = "All maple trees in Michigan Big Woods site")
ggplot(data = bw_growth %>% filter(species == "shrub")) +
geom_sf(aes(col = species), size = 0.5)+
labs(title = "All maple trees in Michigan Big Woods site")
maples <- read_csv("https://rudeboybert.github.io/SDS390/static/data/maples.csv")
Model loos like this: \[ \begin{eqnarray*} y &=& \beta_0 + \beta_{\text{dbh}}\cdot \text{dbh} + \beta_{\text{evergreen_bm}}\cdot \text{evergreen_bm} + \beta_{\text{maple_bm}}\cdot \text{maple_bm} + \beta_{\text{misc_bm}}\cdot \text{misc_bm}\\ && + \beta_{\text{oak}}\cdot \text{oak} + \beta_{\text{short_tree}}\cdot \text{short_tree} + \beta_{\text{shrub}}\cdot \text{shrub} + \epsilon \end{eqnarray*} \]
where
model_maple <- lm(growth ~ dbh + evergreen + maple + misc + oak + short_tree + shrub, data = maples)
model_maple can be written as: \[
\begin{eqnarray*}
y &=& .120157747 + .008636484\cdot \text{dbh} -0.691935999\cdot \text{evergreen_bm} + .079657822\cdot \text{maple_bm} -2.402421803\cdot \text{misc_bm}\\
&& -0.115779874\cdot \text{oak} + .648948573\cdot \text{short_tree} + 2.506065233\cdot \text{shrub} + \epsilon
\end{eqnarray*}
\] Residual analysis Look at model significance
After you’ve selected the model, describe the relationship between \(y\) = growth and \(\text{dbh}\) and \(\text{biomass}\), showing all your work:
model2:\[ \begin{eqnarray*} y &=& .120157747 + .008636484\cdot \text{dbh} -0.691935999\cdot \text{evergreen_bm} + .079657822\cdot \text{maple_bm} -2.402421803\cdot \text{misc_bm}\\ && -0.115779874\cdot \text{oak} + .648948573\cdot \text{short_tree} + 2.506065233\cdot \text{shrub} + \epsilon \end{eqnarray*} \]